library('data.table')
bds <- fread('bds2018.csv')

#bds <- bds[year >= 1985]

pop <- fread('POPTHM.csv', na.strings = '.')
pop$year <- as.numeric(substr(pop$DATE, 1, 4))
pop$DATE <- NULL

bds_e <- fread('bds2018_eage.csv', 
               colClasses = c(year = 'numeric', eage = 'character', firms = 'numeric',estabs = 'numeric',emp = 'numeric'), na.strings = c('(S)'))


bds_e_d <- bds_e[,.(entry_emp = sum(emp*(eage == 'a) 0'))), by = 'year'] 

bds <- merge(bds, bds_e_d, by = 'year')
bds <- merge(bds, pop, by = 'year')

png('Estab_Growth_Per_Capita.png',width = 5, height = 3, units="in", res=300, pointsize=8)
par(lwd = 4)
par(mfrow = c(1,1), family="serif")
par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.3,3.6,1.1,1.1))
plot(bds$year[2:length(bds$year)], 100*diff(log(bds$estabs)) - 100*diff(log(bds$POPTHM)), type = 'l', xlab = 'Year', ylab = 'Log Growth',
     ylim = c(-3.25, 3.25))
rect(c(1980, 1981, 1990, 2001, 2007), par("usr")[3]+.06, c(1981, 1983, 1992, 2002, 2010), par('usr')[4]-.06, col = 'grey', border=NA)
abline(h = 0, lty = 2)
lines(bds$year[2:length(bds$year)], 100*diff(log(bds$estabs)) - 100*diff(log(bds$POPTHM)))
dev.off()

# Entry rates - full sample
png('Establishment_Entry_Time_Series.png',width = 5, height = 3, units="in", res=300, pointsize=8)
par(lwd = 4)
par(mfrow = c(1,1), family="serif")
par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.3,3.6,1.1,1.1))
plot(bds$year, bds$estabs_entry_rate, type = 'l', ylim = c(0, 15), xlab = 'Year', ylab = 'Percent')
rect(c(1980, 1981, 1990, 2001, 2007), par("usr")[3]+.05, c(1981, 1983, 1992, 2002, 2010), par('usr')[4]-.08, col = 'grey', border=NA)
lines(bds$year , bds$entry_emp/bds$emp*100, lty = 2, col = 2)
lines(bds$year, bds$estabs_entry_rate)
legend('bottomleft', lty = c(1,2), col = c(1,2), legend = c('Entry rate','Employment share'))
dev.off()

png('num_entrants_estab_per_cap.png', width = 5, height = 3, units = 'in', res = 300, pointsize = 8)
par(lwd = 1)
par(mfrow = c(1,1), family="serif")
par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.3,3.6,1.1,1.1))
plot(bds$year, bds$estabs_entry/bds$POPTHM, type = 'l',
     xlab = 'Year', ylab = 'Entrants. per 1000 people')
rect(c(1980, 1981, 1990, 2001, 2007), c(0), c(1981, 1983, 1992, 2002, 2010), 1000,
     col = 'grey')
lines(bds$year, bds$estabs_entry/bds$POPTHM, lwd = 4)
dev.off()


png('Estab_Per_Capita.png',width = 5, height = 3, units="in", res=300, pointsize=8)
par(mfrow = c(1,1), family="serif")
par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.3,3.6,1.1,1.1), lwd = 4)
plot(bds$year, bds$estabs/bds$POPTHM, type = 'l', xlim = c(2006, 2018), ylim = c(21, 23),
     xlab = 'Year', ylab = 'Estab. per 1000 people')
dev.off()

par(mfrow = c(1,1), family="serif")
par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.3,3.6,1.1,1.1), lwd = 4)
plot(bds$year, bds$estabs/(bds$estabs[bds$year == 2007])-1, 
     type = 'l',
     xlab = 'Year', ylab = 'Estab. per 1000 people')

#bds_est <- lm(bds$estabs ~ bds$year)
#trend <- lm(bds$estabs ~ bds$year, data = bds_est)


png('Establishment_Entry_And_Exit_Time_Series.png',width = 5, height = 3, units="in", res=300, pointsize=8)
par(lwd = 4)
par(mfrow = c(1,1), family="serif")
par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.3,3.6,1.1,1.1))
plot(bds$year, bds$estabs_entry_rate, type = 'l', ylim = c(7, 17), xlab = 'Year', ylab = 'Percent')
rect(c(1980, 1981, 1990, 2001, 2007), par("usr")[3]+.1, c(1981, 1983, 1992, 2002, 2010), par('usr')[4]-.1, col = 'grey', border=NA)
lines(bds$year , bds$estabs_exit_rate, lty = 2, col = 2)
lines(bds$year, bds$estabs_entry_rate)
legend('bottomleft', lty = c(1,2), col = c(1,2), legend = c('Entry rate','Exit rate'))
dev.off()

# filtered version

png('Establishment_Entry_And_Exit_Time_Series_Filtered.png', width = 5, height = 3, units = 'in', res = 300, pointsize = 8)
par(lwd = 4)
par(mfrow = c(1,1), family="serif")
par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.3,3.6,1.1,1.1))

bds$estabs_entry_rate.trend <- filter(bds$estabs_entry_rate, rep(1/5,5), sides = 1)
bds$estabs_entry_rate.cycle <- bds$estabs_entry_rate - bds$estabs_entry_rate.trend

bds$estabs_exit_rate.trend <- filter(bds$estabs_exit_rate, rep(1/5,5), sides = 1)
bds$estabs_exit_rate.cycle <- bds$estabs_exit_rate - bds$estabs_exit_rate.trend


plot(bds$year, bds$estabs_entry_rate.cycle, type = 'l', ylim = c(-2.5, 2.5), xlab = 'Year', ylab = 'Percent')
rect(c(1980, 1981, 1990, 2001, 2007), par("usr")[3]+.045, c(1981, 1983, 1992, 2002, 2010), par('usr')[4]-.045, col = 'grey', border=NA)
lines(bds$year, bds$estabs_entry_rate.cycle)
lines(bds$year , bds$estabs_exit_rate.cycle, lty = 2, col = 2)
abline(h=0, lty = 2, lwd = 2)
legend('bottomleft', lty = c(1,2), col = c(1,2), legend = c('Entry rate','Exit rate'))
dev.off()

png('Establishments_per_capita_Great_Recession.png',width = 5, height = 3, units="in", res=300, pointsize=8)
par(lwd = 4)
par(mfrow = c(1,1), cex = 1.5, family = 'Arial')
par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.3,3.6,1.1,1.1))
initial_condition <- log(bds$estabs[bds$year == 2007]) - log(bds$POPTHM[bds$year == 2007])
plot(bds$year, 100*(log(bds$estabs) - log(bds$POPTHM) - initial_condition), type = 'l',
     xlim = c(2007, 2018), ylim = c(-8, 0.5), xlab = 'Year', ylab = 'Percent deviation from 2007')
abline(h=0, lty = 2)
dev.off()

png('Employment_per_capita_Great_Recession.png',width = 5, height = 3, units="in", res=300, pointsize=8)
par(lwd = 4)
par(mfrow = c(1,1), family="serif")
par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.3,3.6,1.1,1.1))
initial_condition <- log(bds$emp[bds$year == 2007]) - log(bds$POPTHM[bds$year == 2007])
plot(bds$year, 100*(log(bds$emp) - log(bds$POPTHM) - initial_condition), type = 'l',
     xlim = c(2007, 2018), ylim = c(-8, 0.5), xlab = 'Year', ylab = 'Percent deviation from 2007')
abline(h=0, lty = 2)
dev.off()


png('Establishments_entry_and_exit_per_capita_Great_Recession.png',width = 5, height = 3, units="in", res=300, pointsize=8)
par(lwd = 4)
par(mfrow = c(1,1), family="serif")
par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.3,3.6,1.1,1.1))
initial_condition <- log(bds$estabs_entry[bds$year == 2007]) - log(bds$POPTHM[bds$year == 2007])
plot(bds$year, 100*(log(bds$estabs_entry) - log(bds$POPTHM) - initial_condition), type = 'l',
     xlim = c(2007, 2018), ylim = c(-35, 15), xlab = 'Year', ylab = 'Percent deviation from 2007')
abline(h=0, lty = 2)
#dev.off()

#png('Establishments_exit_per_capita_Great_Recession.png',width = 5, height = 3, units="in", res=300, pointsize=8)
#par(lwd = 4)
#par(mfrow = c(1,1), family="serif")
#par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(3.3,3.6,1.1,1.1))
initial_condition <- log(bds$estabs_exit[bds$year == 2007]) - log(bds$POPTHM[bds$year == 2007])
lines(bds$year, 100*(log(bds$estabs_exit) - log(bds$POPTHM) - initial_condition),
      col = 2, lty = 2, lwd = 4)# type = 'l',
    # xlim = c(2007, 2018), ylim = c(-35, 15), xlab = 'Year', ylab = 'Percent deviation from 2007')
#abline(h=0, lty = 2)
legend('topright', col = 1:2, lty = 1:2, lwd = 4, legend = c("Entry","Exit"))
dev.off()


